
functions {

  real GEVcpdf_lpdf(real x, real mu, real sig, real xsi){
    real z=(x-mu)/sig;
    real tau=exp(-6.0*abs(xsi));
    real lht;		
    if(1+xsi*z<tau){
      lht=-log(tau)/xsi-(z-(tau-1)/xsi)/tau;
    }
    else{
      lht=-log1p(xsi*z)/xsi;
    }
    return (1+xsi)*lht-log(sig);
  }

  real GEVpdf_lpdf(real x, real mu, real sig, real xsi){
    real z=(x-mu)/sig;
    real tau=exp(-6.0*abs(xsi));
    real lht;
    if(1+xsi*z<tau){
      lht=-log(tau)/xsi-(z-(tau-1)/xsi)/tau;
    }
    else{
      lht=-log1p(xsi*z)/xsi;
    }
    return (1+xsi)*lht-log(sig)-exp(lht);
  }

}

data {

    real xi_min;
    real xi_max;
    int<lower=0> T;
    int<lower=0> nobs;
    matrix[nobs,T] y;    
    real sg_xi;
    real sg_alpha;
    real sg_s;
    
}

parameters {

    // Level values for parameters
    real trans_xi_level;
    real ln_s_level;
    real m_level;
    // Innovation standard deviations (unscaled) for random walks
    real<lower=0.0,upper=5.0> g_xi;
    real<lower=0.0,upper=5.0> g_alpha;
    real<lower=0.0,upper=5.0> g_s;
    real<lower=0.0,upper=5.0> g_m;
    // Random walks
    vector [T] h_xi;
    vector [T] h_alpha;
    vector [T] h_s;
    vector [T] h_m;

}

transformed parameters {

    // set sg_m
    real sg_m = exp(ln_s_level);  // This preserves invariance to scale
    // Scale factors for innovations
    real gamma_xi = sg_xi*g_xi/sqrt(T);
    real gamma_alpha = sg_alpha*g_alpha/sqrt(T); 
    real gamma_s = sg_s*g_s/sqrt(T);
    real gamma_m = sg_m*g_m/sqrt(T);
     // vector of zeros
    vector [T] zeros = rep_vector(0.0,T);
    // Random walks
    vector [T] drw_xi = cumulative_sum(gamma_xi*h_xi)
                        -mean(cumulative_sum(gamma_xi*h_xi));
    vector [T] drw_alpha = cumulative_sum(gamma_alpha*h_alpha)
                           -mean(cumulative_sum(gamma_alpha*h_alpha));
    vector [T] drw_s = cumulative_sum(gamma_s*h_s)
                       -mean(cumulative_sum(gamma_s*h_s));
    vector [T] drw_m = cumulative_sum(gamma_m*h_m)
                       -mean(cumulative_sum(gamma_m*h_m));
    //  Construct Basic Parameters
    vector [T] trans_xi = trans_xi_level + drw_xi;  
    vector [T] ln_alpha = drw_alpha;
               // note that level of ln_alpha is fixed at 0
    vector [T] ln_s = ln_s_level + drw_s;
    vector [T] m = m_level + drw_m;
    // Transformed Parameters
    vector [T] alpha = exp(ln_alpha);
    vector [T] s = exp(ln_s);
    // Construct GEV Parameters
    vector [T] xi = xi_min + (xi_max-xi_min)*((exp(trans_xi))
                             ./(1+exp(trans_xi)));
    vector [T] sigma = s.*(alpha.^xi);    
    vector [T] mu = m + (s./xi).*((alpha.^xi)-1);

}

model {

  // Priors
    g_xi ~ exponential(1.0);
    g_alpha ~ exponential(1.0);
    g_s ~ exponential(1.0);
    g_m ~ exponential(1.0);
    
    trans_xi_level ~ normal(0.0,1.5);
    // Flat prior on ln_s_level and m_level (imposed in Stan)

    h_xi ~ std_normal();
    h_alpha ~ std_normal();
    h_s ~ std_normal();
    h_m ~ std_normal();

    for (t in 1:T) {
      for (j in 1:nobs-1) {
        y[j,t] ~ GEVcpdf(mu[t],sigma[t],xi[t]);
      }
      y[nobs,t] ~ GEVpdf(mu[t],sigma[t],xi[t]);
    }
    
}
